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Abstract 

We study an individual based model describing competition in space between two dif- 
ferent alleles. Although the model is similar in spirit to classic models of spatial population 
genetics such as the stepping stone model, here however space is continuous and the total 
density of competing individuals fluctuates due to demographic stochasticity. By means of 
analytics and numerical simulations, we study the behavior of fixation probabilities, fixa- 
tion times, and heterozygosity, in a neutral setting and in cases where the two species can 
compete or cooperate. We argue that this model is a natural starting point to describe more 
complex situations, most notably competition in marine environments where individuals are 
transported by fluid flows. 



1 Introduction 

A mathematical analysis of the fate of mutations in spatially extende d populations has been a 

classic topic of research in population genetics for at lea s t seve nty years (|Fisherl . ll937l : lKolmogorov et al 
19371 : IWrightL Il943t iKimural Il953l : Kimura and Weissl . Il964j ) . This interest has nevertheless in- 
creased recently, as improved sequencing technology allows direct observations of structured 
genetic diversity in space for many different species. 

On the theoretical side, a landmark in this research has bee n the stepping stone model (SSM) 
proposed by Kimura ( Kimura . 19531 : Kimura and Weissl . Il964h . This model considers m islands 
(or "demes"), each having a fixed local population size Ni and arranged along a line or in a 
regular lattice in more than one spatial dimension. The population on each island is made 
up of several species (or alleles) described by, e.g., a Wright-Fisher or Moran process. Spatial 
migration is modeled by assuming that neighboring islands exchange individuals at some given 
rate. 
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It is often convenient to describe the state of the system in terms of the macroscopic density 
of individuals /(x, t) carrying one of the two alleles. In the continuum limit, the macroscopic 
equation governing the time evolution of such density reads 



&/(x,i)=£>V 2 /(x,i) +s/(l-/) 



/(I - /) 



N 



£(x,t) 



(i) 



where N — Ni/a d , a is the lattice spacing between two neighboring islands^, d the spatial 
dimension, and £(x, t) is a Gaussian stochastic process, delta correlated in space and time, 
(£(x, i)£(x',t')) = S(x — x')S(t — t'). Here, / = 1 means an island exclusively populated with 
one allele and / = menas exclusive occupation by the alternative genot ype. The nonlinearit y 
multiplying the noise requires an interpretation in terms of the Ito calculus ( Korolev et all [20091 ). 

However, in many realistic cases, the mechanism of species movemement and range expansion 
is more complicated than a simple diffusion process. For exa mple, recent observ ations on crabs 



colonies along the east coast of north America eastern coast (jPringle et all 120111 ) demonstrated 



how invasion of one allele is controlled by the asymmetrical advection of larvae from north to 
south by a coastal current. The interplay between population genetics and individual movement 
(and transport) can be even more complex in the open ocean, where individu als belonging to 



different planktonic and bacterial species are s tirred and mixed by chaotic flows ([D'Ovidio et al 



20101 IPerlekar et all lioibl iBenzi etail |2012[ ). Of particular interest is the population genetics 
of photosynthctic organisms that control their buoyancy to remain near the surface of an aquatic 
environment. In this case, the advecting flows are effecti vely compressible, leading to popula- 
tion d ensities that overshoot the normal carrying capacity ( Perlekar et al. , 2010l : Pigolotti et al 
2012ft . 



While the SSM can be generalized to include a constant asymmetric diffusion (see i.e. 
(|Pringle et all 120111) 1. ;he extension to more complex fluid environments is quite subtle. One 
of the main underlying assumptions of the SSM - a local population size that does not vary 
either in time nor in space - is quickly violated in acquatic environments where flows create 
inhomogeneities in the total density of individuals. Furthermore, in this case it is less appropri- 
ate to discretize the system in space into demes with a fixed size. In compressible turbulence, 
for example, the dens ity of individuals can be inhomogenous on a wide variety of spatial scales 
([Perlekar et all 12010ft . even inside a single deme (which in the SSM is assumed to be well-mixed). 

In this paper, with the goal of describing population genetics in acquatic environments in 
mind, we introduce a new model in which individuals carrying two different alleles A and B live 
in a continuous space. Their individual densities are allowed to grow and fluctuate, including the 
important possibility of overshooting the natural carrying capacity. Indeed, note that naively 
assuming compressible flows that make / > 1 would lead to an imaginary noise amplitude in Eq . 
([Tt ! The model we study is similar in spirit to the stochastic logistic equation ( Law et all 120031 : 
Hernandez-Garcia and Lopezl . 2003; iBirch and W~aL l200fift . However, in this study we focus 
on competition and cooperation of two species, rather than the stochastic growth of a single 
population. The second differe nce is that previous studies focused on patterns formed by the 
non-local nature of competition ([Hernandez-Garcia and Lopezl . l2003t IBirch and W. AI . I2006T ). In 



x It is convenient to distinguish between jV; (the population inside a single discrete deme of the SSM) and N 
(the corresponding total density of individuals). The former is the quantity used to define the model, while the 
latter determines the amplitude of the noise due to number fluctuations in the continuum formulation of Eq. ((TJ. 
Notice that AT; is a non-dimensional quantity, while N is a density, carrying units of an inverse length to the 
power d in d dimensions. 
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this paper, we mostly focus on the parameter range in which such patterns are not formed and 
a weak noise description in the spirit of Eq. ([1]) is appropriate. 

The phenomenology of such a model, even in the presence of very si mple flows, is very rich due 



to the interplay between population dynamics and fluid advection (see lPigolotti et al.l (|2012l ) for 
some of the consequences in one dimension). For this reason, here we limit ourselves to the case in 
which the flow is absent and individuals move in space in a diffusive way. This simple case allow 
for a systematic comparison with the known results of the SSM. In particular, we show that there 
exists a parameter range where the predictions of our model are consistent with E q. ([1]) and its 
generalization to include competitive exclusion and mutualism ( Korolev and Nelsonl . 12011 ). In 
simple cases, such as when the two species are neutral variants of each other, this correspondence 
can be shown analytically. In more complex cases, the correspondence is explored by means of 
numerical simulations. 

In Sec. [2 we sketch the model of growth, competition and cooperation studied here, which 
leads to the two-species model for allele densities ca(x, t) and c_b(x, t) summarized in Eq. ©. We 
focus on three interesting cases: (1) strictly neutral competitions, (2) a reproductive advantage 
of one species over the other and (3) mutualistic situations where cooperation plays a role. Sec. 
[3] discusses the behavior of our model in the "zero-dimensional" well-mixed case in which the 
population is not structured in space, which allows us to determine limits such that standard 
Wright-Fisher and Moran results for population genetics can be recovered from our more general 
model. Motivated by future applications to population genetics and fluid flows, we then explore 
in Sec. 2] the long-time behavior of our model without fluid advection in one and two spatial 
dimensions. Concluding remarks are presented in Sec. [5] A detailed derivation of our model 
equation is contained in Appendix A. Appendix B shows how conventional stepping stone model 
results can be recovered in certain limits. Appendix C describes a lim it in which a mutualistic 
generalization of the famous Kimura formula for fixation probabilities (jCrow and KimuraL[l97nT) 
is possible. 

2 Model 

Many widely studied models of population genetics in space, the most notable example being 
the stepping stone model, consider individuals carrying different alleles that occupy sites (also 
called "demes") on a lattice. It is commonly assumed that each site is always saturated up to 
its carrying capacity, so that, at each deme, the local population size Ni is constant during the 
dynamics. 

We relax these assumptions by considering discrete individuals Aj carrying different alleles 
(denoted by the index i) and diffusing in continuous space (with a diffusion constant D, for 
simplicity equal for all individuals). Further, we implement population dynamics assuming that 
individuals carrying allele i reproduce at rate jii and die with rates Ay proportional to the 
number of individuals carrying a (possibly) different allele j in a region of spatial size 8 centered 
on their position. For example, in one dimenson (Id), 5 will be an interaction length, while in 
2d it will be an interaction area. In a language borrowed from chemical kinetics, the "reactions" 
we consider are: 



X,- t -4 2Xi (reproduction) 
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Xi + Xj Xi (death by competition) (2) 

In the case of a single species, this s et of reactions is commonly referred to as the birth- 
coagulation process ( Doering et al. , 2003). In this paper, we will focus on the case of two alleles, 



i = A, B. Other reactions could be added to the ones above, for example the possibility that 
an individual can die even in absence of competition, Xi — > 0, or reactions implementing more 
complex biological interactions. We will limit ourself to the biological dynamics embodied in 
(|2|) , which contains minimal ingredients necessary to generate most of the main features present 
in more complicated models. Notice that, in contrast to models such as the Moran process, the 
density of individuals is not fixed but fluctuates both locally and globally. 

In order to make the presentation more compact, we start by discussing the spatially explicit 
version of the model and then discuss the globally well-mixed version as a limiting case. We 
consider the number densities ^(x, t) and ng(x, £), that integrated over a region of space 
yield the (stochastic) number of individuals of species A or B in that region. We will study 
cases in which the number densities are typically large, and consequently define concentrations 
ca(x, i) = n^x, i)/iV and cb(x, t) = Us(x,i)/iV via a constant parameter N, assumed to be of 
the same order of magnitude of ha and ns- This means that, by definition, a constant density 
c = 1 corresponds to a uniform distribution of N individuals in a segment of length 1 in one 
dimension. More generally, in d dimensions, a concentration c(x, t) = 1 will correspond to a total 
number of particles Af — NL d in a system of linear size L. With this choice, the macroscopic 
equations describing the dynamics of the concentrations ca, cb of species A and B read: 



d 2 / c A {fJ-A + Xaaca + Xabcb) c 
—c A = DV c A + c a (^a - Xaaca - Xabcb) + y f 



d 2 i / x , \ . c b {hb + X B aca + Xbbcb) c , (rt , 
—c B = DV c B + c b (^b - Xbaca~ Xbbcb) + y 1, (3) 

where £(x, t) and £'(x, t) are independent Gaussian random variables, delta-correlated in 
space and time, < £ (x, £)£(x, t') >= S(t — t')S(x — x') that should be interpreted according to 
the Ito prescription ( Korolev et al. . 2009). The parameters multiplying the quadratic terms in 



the concentrations are defined in terms of the Ay as \j — NSXij, where S is the interaction 
domain defined above. In the following, we will focus on cases in which the fi^s and the Ay's 
are of the same order of magnitude, so that typical values of the total concentration ca + cb are 
order 1. Under these assumptions, it is useful to note that the quantity 2iV~ 1 = 2(5/ (Ay /Ay ) 
plays here the same role of the genetic diffusion constant in the stepping stone model. In 
particular, 6 is analogous of the lattice spacing, while the denominator on the right hand side 
can be thought as the carrying capacity of each deme. A detailed derivation of Eqs. © is 
presented in Appendix A. If the species densities are well-mixed and we neglect stochastic number 
fluctuations, the deterministic dynamics embodied in Eqs . (131 is a familiar model of growth, 
selection and competition in asexual populations ( Smith! . 119981 ). The four different types of 
dynamcis that emerge depending on the values of the Ay 's are reviewed at the end of this 
section. Our aim here is to understand the rich behaviors possible when both spatial variations 
and number fluctuations are allowed. 

To limit the parameter space, we will consider the following three biologically relevant choices 
for the reaction rates: 
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1. Neutral Theory 

This choice is appropriate when the two biological species (or strains, or mutants and 
wild type alleles) are neutral variants of each other. This means that their growth rates 
and carrying capacities are the same; further, competition with an individual belonging 
to the same species is the same as competition with an individual of the other species. In 
formulas, for Eq. ([3]), a convenient neutral parameter choice is: \la = fJ-B = Aaa = Aab = 
Aba = Abb = M- 

2. Reproductive advantage 

In this setting, we depart from neutrality by allowing for a different reproduction rate of 
species A: fiA = + s) while all the other rates (including the Ay) are equal to fj, as 
in the neutral case. We will study this case to explore the effect of a selective advantage 
of one of the two species on the dynamics of the model. In particular, s > implies a 
selective advantage for A and s < is a disadvantage. Clearly, neutrality is recovered for 
s = 0. 

3. Mutualistic setting 

A simple way to study mutualistic interactions is to assume that the only departure from 
neutrality occurs in the intensity of competition between individuals carrying different 
alleles. In formulas, we have ha = MB = M) Aaa — Abb = M; Aab = m(1 — ca), 
and Aba = m(1 — e b)- The corresponding macroscopic equations are well defined only 
for ca,£b < 1, so that the competition rates Ay are non- negative. We will focus mostly 
on the case ca > and 6r > 0. I n this regime, spatial number fluctuations play an 



important role (jKorolev and Nelsonl . 120111 ) and competition between species is reduced 
(we will interpret this reduction as the effect of mutualistic interactions). Other choice 
could also be of interest, for example ca = and cb < is another way of allowing a 
competitive advantage of A over B (in this case, via enhanced competition rather than 
via a larger reproduction rate). We note finally that ea < 0, £b < corresponds to 
a competitive exclusion model, arising for example when the competing variants secret 
toxins that inhibit the growth ot their competitors. 

In the following, we will measure time in units of a generation time so that /i = 1. A convenient 
choice of the interaction domain is of the order of the average spacing among individuals, 6 — 
1/N, so that Ay = Ay. For simplicity, we will present most of the spatial results for the one- 
dimcnsional version of the model, introducing two-dimensional results only as appropriate. In 
the spatially explicit case, the system is a segment of length L with periodic boundary conditions. 
We will present also two dimensional simulations, where the system is a L x L square, also with 
periodic boundary conditions. 

An even simpler setting we will study to make contact with traditional Moran or Fisher- 
Wright models is the case in which the population can be assumed to be well-mixed, or "zero- 
dimensional" . This limiting case can be easily obtained from the one dimensional case by setting 
5 = L = 1 and ignoring spatial diffusion, since each individual now interacts with every other 
individual in the population. As a consequence of this choice, one now has N = Ay /Ay. In this 
case, the spatial position of the individuals is irrelevant for biological interactions. Clearly, in 
this special case, the individual density is equivalent to the total number of individuals N = J\f. 
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Figure 1: Three illustrative parameters choices in the one dimensional version of the model. In all 
panels D = 1CP 4 and N = 100. The left panel corresponds to the neutral choice in which all rates 
are set to one and initially the two species are randomly distributed with equal concentrations. 
In the center panel, all parameters are set to one except the reproduction rate of allele A (in 
red) which reproduces at a rate (1 + s) with a large selective advantage s = 0.3; in this case, the 
initial fraction of A is 0.1. In the right panel, competition among species is reduced by taking 
ea = £b = 0.7 to enhance mutualism; in this case the two species are randomly distributed with 
equal concentrations in the initial condition. In this case, mutualism insures that the species (or 
alleles) remain spatially inhomogeneous out to very long times. 



In Fig. ([T]), we anticipate some of the results to illustrate the qualitative behaviors that can 
be explored with the three aforementioned parameter choices in one spatial dimension. In the left 
panel, the two alleles are neutral. Despites fluctuation of the total density, the phenomenology is 
similar to that of the Id stepping stone model: as time progresses, the two alleles are demixed and 
fixation occurs by coalescence of the domain boundaries, which can be regarded as annihilating 
random walks. In the central panel, species A (in red) initially constitutes only 10% of the total 
population; however, it has a reproductive advantage over species B. Despite the discreteness 
of individuals and density fluctuations, there are two noisy Fisher waves by which the initial 
minority can take over the entire population. Finally, in the right panel we simulate a case in 
which mixing of the two species is promoted by reducing competition among different alleles. In 
this case, we expect the two species to remain mixed indefinitey in the limit of large system size. 

In the remainder of this section, we introduce some of the concepts we want to investigate 
in the simple case of a well mixed system without number fluc tuations. In deed, intuition about 
mutualistic behavior (and its opposite, competitive exclusion ( Frev . 2010f can be obtained by 
first neglecting both the spatial degrees of free dom and the no i se ter ms in Eq. ([3]). The dynamics 



of this mean field model is then described by (jKorolev et all 120121 ) 



dt 



n A (t) = n A (t) A - \AA n A (t) - \ABn B {t) 



dt 



n B {t) = n B (t) fi B - \ BA n A {t) - ~^ BB n B (t) 



(4) 



Note that the intrinsic carrying capacities (i.e., the steady state densities of one species when 



G 



the other is absent) for this model are 



Na = Ha/^aa, N b = hb/^bb- (5) 

These quantities (we always choose parameters such that Na ~ Nb) play the role of the param- 
eter N that controls stochastic number fluctuations in the derivation of Eq. Q . 

As mentioned above for case 3, an especially interesting situation arises when (1) the two 
species grow at identical rates when the numbers are dilute, so that \ia = Hb = M> (2) also 
the self-competition terms are also identical, Xaa — ^bb] and (3) the effect of cooperation or 
competitive exclusion is then contained exclusively in the cross-interactions 

Aab = XaaO- ~ Aba = Abb(1 — (6) 

With this choice of parameters, the equations for the concentrations ca — tia/Na and cb = 
tib /n b corresponding to system Q read 

d 

— CA = ^CA [1 - CA - CB + £ACb\ 

■^C A = fJ.c B [1-CA-C B + c B c A \ ■ (7) 

The growth rate fi can be absorbed into a redefinition of time; the remaining two parameters 
ca and e b then control the competition under "crowded conditions" , such that the population 
has grown up to satisfy ca + cb ~ 1. If the two variants are nearly identical, differing by, say, 
only a few genes that control secretion of inessential shared amminoacids or a mild toxin specific 
to the other species, it is reasonable to assume je^l, \cb\ "C 1. 

As illustrated in Fig. [2] the deterministic system always has fixed points a t (0, 0), (0, 1) , 
and (1,0). Depending on the parameters, there can also be a fourth fixed point (Smith, 19981) 
located at 

/ * * \ i^Al Ell) / c \ 

(ca> c b) = ; ■ ( 8 ) 

+ Cb — C-AC-B 

When cooperation is favored (ca, £_b > 0, Fig. [5£i) this fixed point is stable, and leads to a steady 
state population fraction /* of A individuals, < /* < 1, with 

/* = -r^-r = — ^— . (9) 

c* A + c* e A + cb 



When competitive exclusion ( Frey . 2010( l is favored {ca,cb < 0, fig. Eb) this fixed point is 



unstable to the attracting fixed points (1, 0) or (0, 1), depending on the initial conditions. Genetic 
demixing, present in strictly neutral systems only due to stochastic number fluctuations, is 
enhanced in this case. Finally, when ca and cb have opposite signs, the fixed point ([8j lies 
outside the biologically relevant domain, and one of the two fixed points (1, 0) or (0, 1) becomes 
globally stable, corresponding to a competitive advantage for one species or the other when the 
population is dense, see fig. [Hi. 

Suppose we now introduce spatial migration and number fluctuations, to recover the full 
model defined by Eq. ([3]). When might we expect fixation probabilities, the global heterozygosity, 
correlation functions etc. to reduce to the familiar results for conventional spatial stepping stone- 
type models with strictly conserved population sizes in every deme? A particularly simple case, 
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Figure 2: Deterministic dynamics of the mutualistic model in zero dimensions without number 
fluctuations. In a), the interactions ca > and es > favor cooperation, and there is a stable 
fixed point {c* A , c* B ) with both densities nonzero. In B), the organisms secrete toxins that impede 
each others growth, so ca < and e_e < and the fixed point (c A , c* B ) is unstable. 



corresponding to the selectively neutral limit tA = (b — 0, is illustrated for a well- mixed system 
in Fig. 3a below: the population grows up and eventually wanders along the line ca + eg = 1, 
until it reaches the absorbing states at (1,0) or (0, 1). A more general situation is ca + £b = 0, 
in which case one variant typically has a simple selective advantage along an invariant subspace 
given by the line ca + cb = 1- If the fluctuations transverse to this line are small (corresponding 
to a large population size), then the usual formulas for fixation probabilities hold, as we show 
later in this paper. In more general situations, however, it is no longer exactly true that the 
population localizes at long times near the straight line ca + eg = 1- Indeed, we have from Eq. 
© that 

c A + c%= * A + €B , (10) 

£A + £b — e^es 

which exceeds 1 along the outwardly bowed incoming trajectories in Fig. and is less than 1 
for the outgoing inwardly curved trajectory in Fig. [2b- However, we do have the approximate 
equality, c* A + c* B « 1, provided \ca + 6b| -C le^esl in Eq. (flQ|) . In this limit, a combination of 
numerical and analytic arguments presented in this paper show that formulas recent l y deri ved for 
mutualistic and competitive exclusion stepping stone models ( Korolev and Nelson . 201 ll ) apply 
to the current model with demographic fluctuations as well, again provided that the overall 
population size TV is sufficiently large. 

What happens if \ia and (jlb are unequal, but €a and e# remain small? In this case, the 
population proportions will certainly change as an initially small population like that in Fig. 3a 
grows to approach the line ca + Cb ~ 1- However, once this line is reached, the subsequent time 
evolution should again be given by stepping stone model results. 
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3 Well-mixed case with number fluctuations 



In this section, we present the results in the simple well- mixed (or "zero-dimensional") version 
of the model. Thus, we keep number fluctuations in Eq. but neglect spatial variations in 
the allele concentrations. 

3.1 Neutral theory 

As previously discussed, it is useful to describe the dynamics of the neutral version of the model 
in the ca vs. cb plane, as depicted in Fig. (0 left). Starting from a dilute initial condition, the 
system evolves rapidly towards to the intrinsic overall carrying capacity given by ca + cb = 1- 
The dynamics is then localized near this line (with fluctuations), until one of the two species 
goes extinct. This behavior contrasts with the Moran process in which the dynamics is rigidly 
confined to the ca + Cb = 1 line, since no fluctuations of the total density are allowed. To 
determine when these fluctuations are small, first note from Eq. ([3]) that in the neutral case the 
total concentration ct = ca + cb obeys a closed equation: 




Figure 3: Neutral dynamics in the well- mixed case, (a) Example of a trajectory in the (ca,cb) 
plane with N — 500. The initial condition is ua = tlb = 20, i.e. a small fraction of a typical 
long time carrying capacity, (b) Decay of the average heterozigosity (H(t)) for different values 
of N. Curves are obtained from simulations of the particle model; each curve is an average over 
10 4 realizations and the error bars are smaller than the size of the lines, (inset) Same curves 
plotted as a function of t/N. Note the data collapse. 



—CT = HC T {l-c T ) + y — £c, (II) 

decoupled from the fraction of species A, f — ca/(ca + cb), where the noise term £ c satisfies 
(£c(t)£c(t')) = S(t — t'). When TV is large, the stationary solution, beside the solution P(c) = 5(c) 
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corresponding to global extinction that will eventually be reacheddon long times of order exp(JV), 
is approximately a Gaussian with average {or) — 1 and variance {c T ) — {ct) 2 = N , which is 
small when N is large. 

We now describe the dynamics of the relative fraction / of individuals carrying allele A, 
f = ca/{ca +cb). The equation for /(f), derived in[B] reads 

I'-M-^TK?*- (12) 

where £/(f) also satisfies (£/(/)£/(£')) — 6{t — f'), and further we have (£/(*)f c (*')) = 0. 
The above equation allows us to analyze the global heterozygosity, which quantifies the loss of 
diversity as time evolves and is defined as the probability H(t) = 2(/(l — /)) that two randomly 
chosen individuals in the population carry different alleles. 

As mentioned above, the equation for cr is independent of / in the neutral case studied here. 
As a result, one can factorize the average over cr and / in the equation for H{t): 



dt 



ff (f) = -^/(i _/)!+£!; 



4</(i-/)><^; 

N ct 



Neglecting the correction of order N~ 2 , we recover for our model with density fluctuations 
the closed equation for H{t) for Fisher- Wright and Moran-type models with a fixed population 
size derived by Kimura, which states that the to tal heterozygosity decays exponentially in well 
mixed neutral systems ([Crow and Kimura , Il970h : 



{H(t)) = H(0) exp(-2f/Af) 
Fig. (J5}d) confirms this exponential behavior in simulations of the model. 



(14) 



3.2 Reproductive advantage 

We now study the probability of one allele to fixate in a well- mixed population of size N 3> 1, 
in the case in which the allele confers a small reproductive advantage s<l. In the same spirit 
as the previous section, we can derive the equation for the relative fraction / = ca/{ca + c b) 
(see[B]). Upon neglecting terms proportional to s/N, the equation in this case reads: 

|/ = M S /d-/) + f/(l-/)^ (15) 

As in Eq. I12[ this result must be supplemented with the equation for the total concentration 
c t — C A + c b ■ Although in the non- neutral case the equation for cr is no longer independent of 
/, one can show that the averages over cr and / factorize up to terms of order s/N or higher 
that can be safely neglected for s <C 1 and AT 3> 1. 

These observations a llows us to recover the formula for the probability of fixation of allele A 
( Crow and KimuraL fl970h . 

2 Notice that, as in the particle model for simplicity death is implemented only via binary reactions (see Eq. 
[2t , the state of global extinction is not accessible in the particle model. Such discrepancy with the macroscopic 
equation could be easily removed by allowing for death even in absence of competition, i.e. the reaction Xi — > 0. 
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1 _ e -sNf 

Pf™ = I _ e -sN ( 16 ) 

where fo is the fraction of individuals carrying allele A once trajectories like that in Fig ([3^) 
reach the line ca + cb — 1 • This result is again similar to Fisher- Wright or Moran models with 
an exactly fixed total population size. Eq. (ITB|) is tested with simulations in Fig. |U 




Figure 4: Fixation in the well-mixed case with reproductive advantage. Probability of fixation 
for different values of s and TV (listed in the figure) for the well-mixed version of the particle 
model, compared with the prediction of Eq. (|16|) . The initial fraction of individuals belonging 
to species A is fo = 0.1, with ct = 1 initially. The inset shows that all curves collapse when 
plotted as a function of sN. These curves are again averages over 10 4 independent realizations. 



3.3 Mutualism 

In the well-mixed limit of the mutualistic model, fixation always occurs at (ca,cb) — (1,0) 
or (0, 1) after a long enough time. However, when the total number of individuals is large, 
this time grows exponentially with N and can easily become inaccessible to experiments (and 
simulations) . As detailed in [C] the quasi-stationary solution where the two cooperating species 
coexist for ca^b > can be seen as a state confined by two potential barriers, one inhibiting 
species A to fixate and the other inhibiting species B to fixate. When N is large, it will be 
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extremely probable that fixation will occur by passing the lowest of these two barriers. In this 
case, an estimate of the time t* needed to reach fixation can be derived by calculating the height 
of the lowest barrier and applying Kramer's escape rate theory. The result is: 



t* 



exp 



Nmm(e 2 A ,e 2 B ) 



(17) 



Figure (JS|) shows a heat map of the total heterozygosity in the (e^, c-b plane for N — 500 after 
5000 generations. The black region is where fixation occured. Green lines are the theoretical 
limits of the apparent coexistence region obtained from Eq. 1171 



<H> 




-0.5 0.5 



Figure 5: Finite-time coexistence in the well-mixed mutualistic model. Average heterozygosity 
in the (6a,£b) plane, with N = 500, in d — dimensions, i.e. for the well-mixed model. 
Simulations are run until a time t — 5000. For each pair of (e^, es) values, after a transient, the 
heterozygosity approaches a quasi-stationary value. The green line limits the region in which 
coexistence up to this time is possible according to the estimate (fTT)l . 



After estimating the fixation time in the mutualistic model, we now ask: what is the fixation 
probability of one of the two alleles? In Appendix C, we show that in the appropriate limit 
the fixation probability for mutualists obeys a formula simi l ar to the result for a stepping stone 
model with fixed total population size (jKorolev and Nelsonl . 120111 ). namely 



«(/o) 



fo iNs(f'-p) 2 



(18) 
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where fo is the initial fraction of allele A. In the limit /* — > oo, s — > 0, with a mutualistic 
effective selective advantage s = f*s fixed, this reduces to the famous Kimura formula discussed 
above 



«(/o) = ; — ivj (19) 



1 _ e -JVS/o 

1 - er 

The formulas above are a good approximation for arbitrary initial conditions only for the 
case of equal initial growth rates fiA = Mb = A*j so that population fractions are approximately 
unchanged prior to reaching the line ca + cb ~ 1- We explore the fixation probabilities in three 
different cases, each having a different definition of selective advantage: 



£a + £b = 0. Unless ea = (b = 0, this corresponds to a selective advantage under crowded 
conditions, such that ca + Cb ~ 1- In the previous section, we discussed how in the 
deterministic limit there are two stable fixed points, {c* A , c* B ) — (1, 0) and (c A , c* B ) = (0, 1), 
while the fixed point with both c* A and c* A nonzero is inaccessible. Fig. [6^ shows the 
fixation probability for c^(i = 0) + cs(t = 0) = 1 and two initial frequencies fo = 0.5, 
fo = 0.1, N = 250, /* = ca/^a + £ s) and effective selective advantage s — (1€a — —\izb- 
The population size N appears through the combination s * N in Eq. 221 so we plot the 
probability versus this rescaled parameter. We obtain excellent agreement between this 
special case of our model and the Kimura formula for the Moran model Eq. (|19l) . 

f-A + e B = s/f*,6A > 0,6b > 0. This corresponds to a mutualistic situation in which 
there is a stable fixed point out in the plane (c* A > 0,c* B > 0). Fig. [BJd shows the fixation 
probability u(fo) versus the initial fraction fo for stochastic Gillespie simulations with 
(A — (b = 0.1 where /* = 0.5 and TV = 100. For comparison, the formula Eq. ([42]) is shown 
as the full drawn line again indicating very good agreement. 

€a + Ef? = —s/f*,eA < 0, cb < 0. This choice corresponds to the competitive exclusion 



(|Frevl . l2010h n which there is an unstable fixed point in the plane (c* A > 0, c* B > 0) and two 



stable fixed points where one of the two species has gone extinct. Fig. [6J; shows Gillespie 
simulations for three cases of €a < 0, eg < and a comparison with the formula Eq. (|42l) 
for the different values of /* (in order to compare this case we take s < in the formula). 

As a further case we consider random initial conditions 71^(0),71b(0) uniformly distributed 
in the square [1,N] x [l,iV], so that the approach to the line ca + cb ~ 1 can play a role as 
well. The initial fraction is now defined as fo = nA(0)/{riA(0) + rts(0)]. Fig. |6ji show the 
corresponding Gillespie simulation results for 200 different initial conditions for three different 
fractions /* = 0.25, 0.5, 0.75. The analytic fixation curves according to Eq. (|4*2"j) are also shown. 
Although the agreement is excellent, we again expect modification when departures from equality 
of the initial growth rates ha and hb are allowed. 



4 One and two dimensions 

Density fluctuations play a more significant role in one and two spatial dimensions, compared 
with the well-mixed situations described in the previous section. For example, depending on 
initial conditions and genetic drift, different alleles can fix in different regions of space; the 
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ultimate fate of the system then depends on how these different regions interact, which in turn 
depends on the choice of the rates. One of the most striking effects of spatial variation in allele 
number and relative proportions is the existence of a regime in which there is a reduction in the 
average carrying capacity, i.e. the average concentration Z is smaller than the value (Z) = 1 
calculated by neglecting fluctuations. The presence of such a regime is illustrated in Fig. ([7J 
as a function of the D and N. Notice that the latter parameter can be properly interpreted 
as an average number of particles per unit length only when N and D are both large enough. 
In the opposite regime, as a consequence of fluctuations, the average number of particles in a 
unit segment is significantly less than N. We quantify this effect by defining an effective average 
carrying capacity (Z) = (n(t))/N where n(t) is the actual number of particles present at time t 
per unit length and the average (. . .) is over time. 

We find significant deviations from the prediction (Z) — 1 when NyDjJi -C 1 . Hcuristically, 
this criterion can be understood as follows. In spatially extended systems, the populations are 
mixed by diffusion. The diffusion scale y/D/n may be considered as an "effective deme size", 
in the sense that individuals within a distance less than y/D/fi are mixed very efficiently over 
a single generation, while individuals separated by a larger distance are spatially decoupled. In 
one dimension, the condition N\J Dj fj, 3> 1 then corresponds to having many individuals in an 
effective deme size. In the opposite limit, this number is small and fluctuations play a much 
more important ro le. This effect is related to the " strong noise limit" of the stochastic Fisher 
equation (see e.g. (|Doering et all 120031 iBerti et all koOTt iHalfatschek and Korolevlliooi )). We 
remark that in this regime, the assumptions needed to derive Eq.Q from the particle model arc 
violated and significant deviations between the particle simulations and the macroscopic theory 
are expected. For this reason, we will restrict our analysis here to the "weak noise" case in which 



4.1 Neutral theory 

To study how fixation occurs in space, we now discuss the behavior of the spatial heterozigosity 
H(x, t) defined as the probability of two individuals at distance x and time t to carry different 
alleles. In the neutral stepping stone model with a fixed population size in each deme, H (x, t) 
obeys a closed equation: 



d t H(x,t) = DV 2 H-^H6(x) 
In one dimension, such equation can be solved explicitly: 



H(x,t)=H 



2 

1 

N 



dt'- 



erf 



o y/8nD(t - t>) 



' 8D(t-t') 1N' 2 D 



(20) 



(21) 



where Hq is the initial heterozigosity, equal to one half if the two variants are well mixed and 
equally populated at time t = 0. Eg . (12H can be de rived directly from the stochastic Fisher 
equation (JTJ with s = (see, e.g., Korolev et al. ( 20091 )). 

We define the heterozygosity in our off-lattice particle simulations with growth and compe- 
tition from the statistics of intcrparticle distances. In particular, at a given time t, we compute 
all distances between pairs of individals. Upon introducing a bin size h, the function H(r, t) is 
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then defined as the ratio between the number of pairs carrying different alleles at a separation 
between r and r + h, divided by the total number of pairs of all types in the same range of 
separation. For simplicity, we always took the bin size h equal to the interaction distance 5. 

In the limit Ny/D/n 3> 1, the spatial heterozigosity obtained by simulations of the neutral 
off-lattice model shows a remarkable agreement with Eq. (|21|. as shown in Fig. (jHJ. This 
corrrespondence arises because the relative fraction of allele A, f(x,t) = ca/{ca + cb), obeys 
a very similar equations as discussed in the mean field case. In Appendix B, we show that the 
only effect of density fluctuation is an additional effective advection term in th e equation for djj , 



equal to 2D(V log ct) ■ V/. The appearance of such term was already found in lVlad et al.l ([20041 ) 
in a deterministic version of the model described here. In our case, one can show that since Ct 
obeys a decoupled equation in the neutral case, such term will not affect the equation for the 
heterozygosity. Indeed, numerical simulation shows that the average spatial heterozigosity in 
the model reproduces that of the stepping stone model even in the limit of very high diffusivity 
shown in Fig. [51 panel b. Panels (c) shows that similar agreements arise comparing numerical 
integration of Eq. (|20p with our off-lattice simulations in two dimensions. At variance with 
one dimension, where the local heterozygosity H(0,i) decays at long times as i -1 / 2 , in two 
dimension the decay is much slower, H(0,t) l/ln(i). Such slow logarithmic decay is confirmed 
in simulations in panel (d). 

4.2 Reproductive advantage 



In one spatial dimension, an analogue of Kimura's formula (|16|) ([Crow and Kimura , 1970) for 



the fix ation probability has been derived from the stochastic Fisher equation bv lDoering et al 



(2003) 



Pfix = 1 - exp 



sN / dx f(x,t = 0) 



(22) 



where f(x, t = 0) is the initial spatial distribution of the fraction of species A. Remarkably, 
the one dimensional fixation probability is independent of the spatial diffusion constant. We 
tested this prediction in Fig. (JSk), left panel, for our model when species A enjoys a repro- 
ductive advantage s. There are again no appreciable differences between the simulations of our 
more general growth model and the theoretical prediction for the stepping stone model, over a 
wide range of diffusion constants. This agreement is expected, given the approximate mapping 
onto a stepping ston e model embodied in Eq. (B.3) of Appendix B. While the result (|22|) by 
iDoering et aD (|2003h was derived in one dimension, we conjectured that the same formula holds 
in two dimensions. Indeed, a straighforward generalization of Eq. (f22|) predicts well the fixation 
probability in two dimensions, as shown in simulations in Fig. right panel. 

4.3 Mutualism 

We now set \ia — — but allow variable intersp ecific competition coefficients can vary 



in one and two dimensions. iKorolev and Nelson! ([20111 ) recently demonstrated how for a mu- 



tualistic stepping stone model with fixed deme size in one dimension, there is a region in the 
(cAjCb) parameter space in which (in limit of an infinite system size, L — > oo) fixation never 
occurs, as sketched in Fig. [TOl panel a). This behavior differs dramatically from the well-mixed 
zero dimensional case, for wich fixation is inevitable, with a fixation time t* (eA,es, N) given 
approximately by Eq. (fl7|) . 
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We fix parameters as \i = 1, D = 0.02 and N = 30. To explore the behavior of our model, 
we performed simulations along the paths shown as dashed lines in panel a) of Fig. (fTU| . Panel 
b) shows the time evolution of the local heterozigosity H (0, t) along the line ea = eb- For small 
values of £/i = £s > 0, the heterozigosity decays in a similar fashion (roughly as 1/Vt) as in the 
neutral case ea = (b = 0. For higher values, the local heterozygosity eventually levels off at a 
nonzero value, implying that fixation will never occur. 

The presence of a mutualistic regime where the system remains mixed forever is even more 
evident in Fig. [TfJl panel c), where we plot along the cuts at constant ea + eb the long-time 
average of the fraction of the first allele (/) as a function of the difference ea — eb ■ Along the cuts 
that do not cross the mutualistic region, (/) is either or 1 as one of the two alleles always fixes. 
A special case arises for ea = eb , where each of the two alleles has a chance of being fixated equal 
to its relative abundance in the initial condition, so that (/) = Jo- Conversely, (/) has a non- 
trivial behavior along the line ca + £s = 0.4. Upon varying the parameter p = ea — cb, we find 
a wh ole range of values in which fixation does not occur. As discussed in ( Korolev and Nelsonl . 



20111 ). the two lines of critical points shown in (a) are in the directed percolation universality 



class. The behavior of the density close to this critical point is described by a universal exponent, 
ca ~ (eA — £ c )^ f where the e xpected exponent is (3 ~ 0.2765 and e c is the value of ea at the critical 



point (see e.g. lOdorl (|2004j)). Fig. [TU1 panel d) confirms the power law behavior close to one of 



the critical points on the cut ea + eb = 0.4. Finally, in panels e) and f) we show simulations on 
the two dimensional mutualistic model. Mutualism in 2d is computationally challenging and, to 
the best of our knowledge, has not been studied systematically in the literature. Although we 
did not obtain the full phase diagram, our simulations suggest a scenario similar to the Id case. 
In particular, the heterozygosity H(x — 0,i) along the cut ea — £b displays a transition from 
a regime in which it seems to decay logarithmically (as in the 2d neutral version of the model) 
to a regime in which fixation does not occur. Furthermore, the cut at ea + eb = 0.4 shown in 
panel f) reveals a directed-percolation-like transition, qualitatively similar to that in panel c). 



5 Conclusions 



The understanding of growth, competition, cooperation and diffusion in space in individual 

based models has b een subje ct of intense study, in contexts as dive rse as popul ation genetics 

dBarton et al. . 20021). ecology dLaw et all 2003 ; Birch and W. A . 2006f ) and physics ( Hernandez-Garcia and Lopezl . 



2003; iBerti et all l2007t iKorolev et all 120091) . A main focus has been to explore the regime in 



which discreteness effects are such that individual based simulations differ significatively from 
the behavior of macroscopic continuum equations, such as the Fisher equation or its stochastic 
variant. 

In this paper, we have explored competition and cooperation between two different alleles 
when the total population size is not constrained. We have deliberately focused on the weak 
noise limit by choosing carrying capacities and diffusion constants such that there is a good 
agreement between the outcome of the macroscopic Langevin equations and the individual based 
simulations. 

We have shown that, in certain limits, one can draw an explicit correspondence with stepping 
stone- like models in which the total density of individuals is kept fixed at every deme, by studying 
the relative fraction of one of the two species. In the neutral case, the fluctuating total density 
appear in the equation for the relative fraction, but its fluctuations average out in the equation 
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for mean quantities such as the heterozygosity. The correspondence between stepping stone 
models and our generalized off-lattice model with additional fluctuations in the overall density 
was confirmed by individual based simulations. In non-neutral settings, the total density does 
not obey a closed equation and such exact correspondence can not be drawn. However, we have 
shown how, when the departure from neutrality is not severe (small s or small ca and e B ), the 
corrections due to density fluctuations can be safely neglected and the predictions of constant- 
density models are still reproduced with accuracy. The issue we address here is a more subtle 
dynamical version of justifying the neglect of number fluctuations in the grand canonical ensemble 
as compared to the canonical ensemble in equilibrium statistical mechanics. We conclude that 
the model we present here is a natural candidate to study situations in which the total density 
of individuals can vary greatly from the back ground carrying capa city due to external forces, 



such as turbulence or compressible fluid flows (|Pigolotti et all 120121) 
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A Derivation of the macroscopic equations 



In this section we present an explicit derivation of the coupled stochastic macroscopic equations 
for Ca{x, t) and c B (x,t), Eq.©, from the microscopic rate reactions ©. The formalism we wi ll 
follow is that of the chemical master equation, as presented for example by iGardiner ( 2004 ). 
which in turn may be conside red as a spatial generalization of the Kramers-Moyal expansion 
( Gardinerl . l200i iRiskenl . Il989t) . 

As discussed in the section [5J we consider interacting individuals in a volume equal to L d in d 
dimensions. In particular, competition occurs when individuals are within a small volume 5 (for 
details on the implementation of the individual-based dynamics see iPerlekar et al. (2011)). We 
can then discretize the system in cells of size 5 and start the derivation from the master equation 
governing the time evolution of the probability the numbers of particles {nf, n^} of type A and 
B in each cell, labeled by the index j. We first define as R^(±l,n^,n^) and W B (±l 7 nf,n^) 
the rates at which the populations of type A (or B) increase/decrease by one individual in a 
specific box, given that the numbers are currently nf and . The expression for these rates 
are then: 



W A (+l,nf,nf) = fj, A nf 

W A (-1, nf, nf) = X AA nf(nf - 1) + X AB nfnf 
W B {+l,nf,nf) = n A nf 

W B (-l,nf,nf) = X BA nfnf + X BB nf(nf - 1). (23) 
The master equation governing the evolution of the full probability distribution P({nf, nf}, t) 
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for all possible box occupation numbers {nf,nf} then reads: 



3 > 3 

d 
dt 
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diffusion terms, 














(24) 



where the diffusion terms allow for the stochastic exchange of particles between neighboring 
boxes. Although we did not write them explicitly, they reduce to discrete approximations to 
Laplace operator. Indeed, we replace them with Laplacians in the continuous space limit at the 
end of the calculation. 

The next step in the derivation is to perform a Kramers- Moyal expansion ( Risken . 19891) in 
each of the boxes, which leads to 



(25) 



j k=l 



with 



A,B , A B\ ii A,B / a A.B\kjxrA,B i \ A,B -A B\ 

a k { n f,nf)= / dAn j K An j ) W ( A? V >n,J iTlj), 



(26) 



and where the integral over An accounts for the possible jump processes (+1 and —1 in our case). 
Finally, truncating the Kramers-Moyal expansion up to second order in the derivatives leads to 
a Fokker-Planck equation for P{n^,nf}. It is convenient to write directly the equivalent but 
somewhat simpler system of Langevin equations corresponding to this Fokker-Planck description, 
namely: 

dn{ 



dt 
dn B 



rij (fi A - ^AArij - XabtLj ) + diffusion + y n A (fi A + ^AAnf + X A Bnf)£ 



dt 



— nf(fi B - ^BA^f - A BB nf ) + diffusion 



nf(fi B + \ BA nf + \ BB nf)if. 



(27) 



In the above system of equations, the £'s are delta-correlated unit variance Gaussian pro- 
cesses, < £,f {t)£P {t r ) >— SjiS km 5(t — t'). The multiplicat i ve noise in the equa tion must be 
interpreted according to the Ito prescription ( Gardiner , 2004 ; Korolev et al. . 20091) . In principle, 
the diffusion terms in (|2~4"|) would contribute to the noise term. However, on e can show that this 
contribution can be neglected if the size of the cells is sufficiently large (see Gardiner ( 20041 )). 

From Eqs. (l27l) one can finally derive Eqs. (|3]) by: 

1. Taking (formally) the limit 5 — > 0. In such a way the number densities of individuals are 
continuous functions of the coordinate x, n^(x, t) and ns(x, t). 
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2. Defining rescaled, macroscopic rates of binary reactions, 



Ay = NSXij (28) 

3. Defining the macroscopic concentrations of individuals ca,b{^, t) — t)/N. 

The convenience of introducing the macroscopic binary reaction rates Ay in step (2) is that 
the microscopic interaction radius S does not appear in the macroscopic system of equations ([3]). 
At the same time, we introduced a parameter N = Ay/((SAy) that, as clear from step (3) in the 
above procedure, sets the typical number density of particles corresponding to a macroscopic 
concentration c(x, t) = 1. Such parameter does not appear in the deterministic drift terms of 
the equation but only in the noise terms, whose amplitude vanishes for TV — > oo. It is worthwhile 
remarking that, while we followed here the Kramers-Moyal expansion procedure, in the Van 
Kampen formalism the parameter TV -1 is the relevant expansion parameter which is assumed to 
be small (|Riskenl . Il989t bardineii |2004 . 



We remark that through the paper we presented only results of the particle models, corre- 
sponding to given parameter choices in the macroscopic equations © . Equation (|2"51) can be seen 
as defining the mapping between the parameters used in the particle simulations (the interaction 
domain 5 and the microscopic binary rates Ay's) and those appearing in the macroscopic de- 
scription (TV and the set of Ay's) The same relation can be used for the reverse task, i.e. finding 
microscopic parameters 5 and Ay's corresponding to given N and Ay's. Clearly this mapping is 
not univoquely determined, but has one degree of freedom. As sketched in Sec. @, we fixed this 
degree of freedom in two different ways in the well-mixed version of the model and in the d > 0, 
spatially explicit simulations. In particular, in d = we chose 8 = 1, so that the microscopic 
binary reaction rates are N times smaller than the macroscopic ones, Ay = Ay/iV. In this case, 
it is crucial to set the system size L = 1 so that all particles interact with all other particles. 
Instead, in d > we chose the interaction domain S = 1/N, so that the microscopic and macro- 
scopic re action rates ar e ident ical, Ay = Ay. Further details on the simulation schemes can be 
found in lPerlekar et al. l l|201ll) . 

We conclude this Appendix by noting that the continuous space limit is a formal one, and 
cannot be performed in a rigorous way. One of several subtleties is that neglect of the diffusive 
contribution to the noise variance requires a finite value of 5, so that the limit of vanishingly 
small interaction range cannot be taken in a strict sense. Thus, Eq. (j2~3| should be regarded 
as a continuum shorthand notation: In practice, we always simulate equations such as Eq. (|27p 
on a lattice of finite size, and require a smoothly varying total density of particles. When this 
assumption is invalid, the macroscopic description can break down, as briefly discussed in the 
beginning of section (j4]) for the problem of the reduction in the total number of particles for 
d = 1 and d = 2. 



B Appendix: equations for the relative fraction of one 
species 

The correspondence between the growth model presented here and the stepping stone model with 
Fisher- Wright or M oran dynamics, or the equivalent stoch astic Fisher-Kolmogorov-Petrovsky- 
Piscounov equation ( Fisher . 19371 iKolmogorov et al. . 1937 ) can be illuminated by constructing 
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the dynamical equation for the relative fraction of species A, f = ca/ct with ct = ca + c B . A 
dynamical equation for / can be derived with help of the Ito calculus: upon writing the system 
of equation (j3J as: 

a A {cA,c B ) + a A (cA,c B )^(^t) 

a B (c A ,c B ) + (J B {c A ,c B )£,'(x,t) (29) 

where the diffusive Laplacian terms are included into a a , ol b . The equation for the ^-fraction 
/ then reads 

j t f = a A d A f + a B d B f + ^a 2 A (d A f) 2 + <7 2 B (d B f) 2 Z + 

+ ^-d AA f + ^fd BB f, (30) 

where we used the abbreviated notation 8a = d CA , Baa = 9% and so on. Inserting the complete 
set of equations [3] into ([3D]) leads to a lengthy expression for the dynamics of/. However, with the 
choice of parameters we made to discuss a reproductive advantage (this reduces to the neutral 
case for s — 0), Eq. (l3D|) simplifies to 

d 

-g- t ! = DV 2 / + 2£>V(logc T )-V/ + 

Upon neglecting small contributions of order s/N < 1 in the last two terms and neglecting 
fluctuations in the total density (i.e. imposing ct — 1), we recover exactly the equation ([1]) 
governing the macroscopic dynamics of the stepping stone model. 
Repeating the calculation in the case of mutualism yields: 



d / s 



If 
dt J 



DV 2 f + 2£V(logc T ) ■ V/ + M /(l - f)[e A - (e A + e B )f] + 
/*/(/-!), 



N 



[e A (f-l)+e B f} + \ 



M/(W) [(^-)-e A {l-f?-e B P 



N 



< (32) 



Upon neglecting, similar to the case of reproductive advantage, terms order ca, B /N, and again 
neglecting fluctuations away from the line cr(x, t) = ca(x, t) + c B (x.,t) = 1, we recover th e 
continuum limit of the mutualistic stepping stone model treated bv lKorolev and Nelson! (|2011l ). 
namely 

|/ = DW 2 f + s f(l /)(/* - /) + \j ^ f{l N ~ /} e(x, t), (33) 
where s = fi(e A + e B ) and /* = e A /(eA + es)- 
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C Appendix: Fixation times for the mutualistic model in 
the well-mixed case 



To estimate the average fixation time for the mutualistic model in the well- mixed limit, we start 
from Eq. [35] Upon neglecting terms order €a/N, cb/N and also neglecting density fluctuations 
by imposing ct = ca + cb = 1, we obtain: 



dt J 



fif(l-f)[e A -(e A + e B )f} + 



2/*/(l - /) 



N 



(34) 



The dynamics of such equation will reach one of the two absorbing states at / = or / = 1 for 
long enough times. However, these times can be very long when N is large: a time-independent 
metastable probability distributi on exists before the absorbing states are reached, which can be 
written using potential methods ( Gardineil [2004') as 



V{f) oc e- y{f) 



where the potential V is given by 



V(f) = -Nf 



{c A + e B ) 



M/(i-/)] 



(35) 



(36) 



where the first term is analogous to a potential energy and the second resembles an entropy. In 
the large N limit, the potential has a minimum at f c w ca/{za + €b) and two maxima, one at 
/~ sa 1/(Nca) and one at / + « 1 - 1/(Ncb)- By evaluating the potential at these points one 
can estimate the lifetime of the metastable state from the height of the two potential barriers. 
To the leading order in N, the smallest barrier is given by: 



AV 



Nmm(e 2 A ,e%) 



(37) 



Finally, we assume that fixation always occurs via the smallest barrier. With this assumption, 
the time needed to escape the potential minimum to one of the absorbing state can be simply 
estimated from Kramer's escape rate theory as t* ~ exp(AV), which leads to Eq. (fP7|) . 

We now discuss the fixation probability in zero dimensions. The Kolmogorov backward 
equation corresponding to the stochastic differential equation (|34|) . when interpreted using the 
Ito calculus, reads: 



du(p,t) 
dt 



1 d 2 



§p(i-p)(r 



. d 



(38) 



where u(p,t) is the probability that species A has fixed at time t > given that it was present 
with frequency p at time t = 0. We have set /* = ^/(e^ + es), and defined the mutualistic 
advantage s = /j,(ca + £s)- 

Note that Eq. includes the original Kimura problem of two non-interacting species as a 
special case, in the limit /* — > oo, s — > with the selective advantage given by the fixed product, 
s = f*s = fie a- We now define the long time fixation probability for the initial condition p = /q 
as 

limu(/o,t) = u(/ ) (39) 
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Upon assuming a steady state arises at long times, we have from Eq. 

^u'( P ) = Ns(f*-p)u'(p) (40) 

which leads to 

u'{p) = C e ^(/*-p) 2 (41) 

With boun dary conditions u(0) = 0, u(l) — 1, we integrate once more to obtain the fixation 
probability (jKorolev and Nelson , l201lh 

u(f ) = ^Vh , (42) 

a closed form expression in terms of the parameters /oi f*, N and s. It is straightforward to show 
that in the limit /* — > 0, s — > with s = f*s fixed (two noninteracting species with a selective 
advantage s) we recover Kimura's famous formula for the fixation probability, Eq. (1191) . 
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Figure 6: Fixation probabilities u(s, fo, -/V)in the mutualistic model. Full curves show the ana- 
lytical results from Eq. (|42"T) with initial fraction fo = ua (0)/(7M.(0) + ?ib(0). a) Competitive 
exclusion: Simulations with + e B — with —0.08 < 6a < 0.2, N = 250: Red curve: fo = 0.1, 
blue curve f = 0.5. Green and purple curves: Eq. HU) with -0.08 < s < 0.2, N = 250, 
fo = 0.1 (green), fo = 0.5 (purple). The curves are plotted against the scaling variable 
s * N for different initial frequencies fo- Here (and also in b),c)) the initial condition is 
chosen on the line n^(0) + n B (0) = N. b) Mutualism: Green curve: simulations with 
(A = e B = 0.1, N = 100, /* = 0.5. The fixation probability u is plotted versus the initial fraction 
fo- Red curve: Fixation formula p^|) with N = 100, s = 0.1,/* = 0.5. c) Coordination game 
with an unstable fixed point /*: Green; purple; orange curves: simulations with = —0.05, cb — 
-0.15(/* = 0.25); 6a = -0.10, e B = -0.10(/* = 0.50); e A = -0.15, e B = -0.05(/* = 0.75) . 
Red; blue; cyan curves: Fixation formula ® with iV = 100,s//* = -0.2,/* = 0.25; 0.5; 0.75. 
d) Mutualism with stochastic initial conditions. Simulations with initial conditions u^i(0), tib(0) 
are uniformly distributed in the plane of size N x N with N — 100. For each random initial 
condition, which fixes the value of /o, the fixation probability is averaged over 500 independent 
Gillespie simulations resulting in u(fo). Cyan points: ca = 0.05, €b = 0.15,/* = 0.25; blue 
points: e A = 0.10, en = 0.10,/* = 0.5; red points: e A = 0.15, e B = 0.05/* = 0.75. Full curves: 
fixation formula (02} with N = 100, s/f* = 0.2: Brown: /* = 0.25; purple: /* = 0.25; green 
/* = 0.75. 
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D 

Figure 7: Reduction of the carrying capacity in the neutral model in Id and 2c?. (a) Reduction 
of the total carrying capacity Z = (c-a + cb) hi the (-D, N )-plane. The system is one dimensional 
and we adopted the neutral choice of parameters (see Section 2). The blue line is the theoretical 
condition NyD/fj, = 1 . (b) Comparison of the carrying capacity reduction in Id and 2d, as a 
function of the nondimensional parameter N 1 / d y/~D j ' [i where d is the spatial dimension. 



25 



D = 0.00001, N = 500 2D heterozigosity 




Figure 8: Heterozigosity in the Id and 2d neutral case. Behavior of heterozygosity correlation 
function for the neutral off-lattice model of growth and competition, (a) ID simulations at low 
diffusivity, D = 10~ 5 and (b) high diffusivity, D — 0.1. In the top case, the system size is L = 1 
while in the bottom case the system size is L — 100. In both cases we find excellent agreement 
with the prediction of formula (|21[) . (c) Neutral heterozigosity in 2d, compared with a numerical 
integration of Eq. (|20|) . (d) Behavior of the local heterozigosity H (x = 0, t) as a function of time 
in 2D, showing the logarithmic decay H(x — 0,t) ~ l/\n(t). 
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Figure 9: Probabilities of fixation in the presence of a reproductive advantage. The two panels 
show (a) one spatial dimension and (b) two dimensions, as a function of the selective advantage 
s, for different va lues of the diffusion c onstant D. Our Id results are compared the results with 
the prediction of iDoering et al.l ( 2003I ). In Id, parameters are N — 500 and the initial fraction 
of species A is fo = 0.01, randomly distributed on the unit interval. The 2d simulations were 
conducted on a square domain of unit area and the parameters TV = 16384 and the initial fraction 
of species A is fo = 0.01 were kept fixed. The solid line is our conjectured generalization of Eq. 
(l22l) to two dimensions. 
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Figure 10: Mutualism in Id and 2d. a) Phase diagram of the mutualistic model in Id. The 
mutualistic region, where global fixation never occurs in an infinite system, is colored in red. 
Dashed lines denote the cuts relevant to data in the other panels, b) Behavior of the local 
heterozigosity H(Q,t) along the cut ca = es- A nonzero long time asymptote implies that 
fixation never occurs, c) Average concentration of allele A, < ca >, along three cuts such that 
£a + £s = const.. When tA + £b is sufficiently large and positive, (/) varies smoothly between 
and 1 when traversing the red region in (a). For both €a + £s = and £a + £b negative, there 
is an abrupt jump in (/) from to 1 when ea = £_b- In this sense, the dashed diagonal line 
below the cusp in (a) is like a first order phase transition. In all figures, parameters are: [1=1, 
D = 0.02, N = 30 and L — 2000 so that on average there are 6 • 10 4 individuals in the system, d) 
Logarithmic plot of the density of A close to the critical point. A power law with the expected 
directed percolation exponent, f(x) oc x 13 , /? » 0.2765 is shown for comparison, e) Behavior of 
the local heterozigosity H(0,t) in 2d along the,jhne ea — ^b- A phenomenology similar to the 
Id case is observed, f) Transition along the diagonal cut ea + es = 0.4 in 2d, again showing a 
similar behavior to the Id case. 



